#include <ctype.h>
#include <math.h>
#include <errno.h>
#include <string.h>
#include "dosearch.h"
#include "parse.h"

#define INDEXFILENUM 6
#define GETSTRING(x) (fscanf(readfile,"%s",x))
#define GETSTRING2(x) (fscanf(readfile,"%s",x),x)
#define GETLINE(x) (fscanf(readfile,"%255[^\n]",x))
#define GOTONUM (fscanf(readfile,"%*[^1234567890]"))
#define GETCURDIR (strcmp(dirName,"") ? dirName : ".")
#define FSCANLINE(x) fscanf(readfile,"%1023[^\n]\n",x)

#define NO_FORK

extern void reset_bss_results(void);
extern const gchar* new_bss_subdir(void);

static void prepare_bss_results_dir(void);

extern char* doctor_bes_name(char* name);
void blat_search(struct bss_instance* pbsi, FPC_TREE* pqtree, FPC_TREE* pdbtree);
void blast_search(struct bss_instance* pbsi, FPC_TREE* pqtree, FPC_TREE* pdbtree);
extern int Zbatch_flag;

extern GtkWidget* cutoff_entry;
extern GtkWidget* additional_entry;

extern char batch_bss_file[MAXPATHLEN];

int g_gb_msg_shown;
int g_seqctg_msg_shown;
int g_length_warn;
int g_remove_bss_msg_shown;
int g_match_msg_count;

void fpc_exec(char* cmd)
{
    printf("execute: %s\n",cmd);fflush(0);
    system(cmd);
}


void collect_filenames(char* dir, char* suffix, GArray* out)
{
   char dir_str[MAXPATHLEN];
   char fil_str[MAXPATHLEN];
   char *ptr;
   struct dirent *dp;
   DIR *dirp;
   struct stat tmptr;
   FILE *readfile;
   char aString[256];
   gchar* tempstr;

   strcpy(dir_str,dir);
   if((dirp=opendir(dir_str))==NULL){
      fprintf(stderr,"Unable to open %s\n",dir_str);fflush(0);
      return;
   }
   while((dp=readdir(dirp))!=NULL){
      if(strcmp(dp->d_name,".") && strcmp(dp->d_name,"..")){
         sprintf(fil_str,"%s/%s",dir_str,dp->d_name);
         stat(fil_str,&tmptr);
         if((tmptr.st_mode & S_IFMT)==S_IFREG){
            ptr=rindex(dp->d_name,'.');
            if((ptr==NULL) && suffix[0]!='\0') continue;
            else{
              ptr++;
              if((suffix[0]!='\0') && strcmp(ptr, suffix)) continue;
            }

	    if((readfile=fopen(fil_str,"r"))!=NULL){
	       GETLINE(aString);
	       if(strstr(aString,">")!=NULL){
                tempstr = g_strdup(dp->d_name);
                g_array_append_val(out,tempstr);
	       }
	       else printf("%s is formatted incorrectly.  Ignored.\n",fil_str);
	    }
	    else{
	       fprintf(stderr,"Cannot open %s: ",fil_str);
	       perror("FPC error");
	       return;
	    }
	    fclose(readfile);
         }
      }
   }
   closedir(dirp);

}


void delete_seq_data(void** ppdata)
{
    assert(ppdata);
    assert (*ppdata);
    free(*ppdata);
    *ppdata = 0;
}


/* Match the query to a clone as follows:
*
* Query <prefix><clone><suffix> matches <clone>
* Query <clone> matches <clone><suffix>
*
* In the first case, the prefix is stored as an optimization.
* (Matching with no prefix is already fast using binary search)
*/

int one_is_contained(char* p1, char* p2)
{
    if (strlen(p1) >= strlen(p2) && 
            0 == strncmp(p1,p2,strlen(p2)))
    {
        return 1 ;
    }
    else if (0 == strncmp(p2,p1,strlen(p1)))
    {
        return 1;
    }
    return 0;
} 

int match_query_to_clone(char* _query, int* pcidx, struct bss_instance* pbsi)
{
    CLONE* clp;
    int result = 0;
    char* query = strdup(_query);
    char* ptr1 = query;
    int i;
    int c;
    int minlen = 2;
   
    if (pbsi->matched_clone == 1 )
    {
        /* we matched one before and we now hope the same prefix will work */
        if (pbsi->prefix_len <= strlen(query))
        {
            ptr1 += pbsi->prefix_len;
            //ptr2 = ptr1 + strlen(ptr1) - pbsi->suffix_len;
            // *ptr2 = 0;
            if (fppFind(acedata, ptr1, &c, cloneOrder))
            {
                *pcidx = c;
                result = 1;
            } 
            else if (c >= 0) 
            {
                /* Because of the way the binary search works, either this
                     or a neighbor  is the best match to the query. If the 
                    query is totally contained in the clone, or vice versa,
                     we will use the match. */

                for (i = -1; i <= 1; i++)
                {
                    if (c+i < 0) continue;
                    if (c+i >= arrayMax(acedata)) continue;
                    clp = arrp(acedata,c+i,CLONE);
                    if (one_is_contained(clp->clone,ptr1))
                    {
                        result = 1;
                        *pcidx = c+i;
                        break;
                    }
                }
            }  
        }           
    }

    if (result == 0)
    {      
        /* Peel off the prefixes, looking for a match. 
            Use the binary search to speed it up. 
            Store the prefix if found. */
        for (ptr1 = query; ptr1 < query + strlen(query) && strlen(ptr1) >= minlen; ptr1++)
        {
            fppFind(acedata, ptr1, &c, cloneOrder);
            if (c >= 0)
            {
                for (i = -1; i <= 1; i++)
                {
                    if (c+i < 0) continue;
                    if (c+i >= arrayMax(acedata)) continue;
                    clp = arrp(acedata,c+i,CLONE);
                    if (one_is_contained(clp->clone,ptr1))
                    {
                        result = 1;
                        *pcidx = c+i;
                        pbsi->matched_clone = 1;
                        pbsi->prefix_len = ptr1 - query;
                        if (g_match_msg_count <= 3)
                        {
                            printf("Matched target sequence %s to clone %s\n",query,clp->clone);fflush(0);
                            printf("Prefix length:%d\n",pbsi->prefix_len);fflush(0);
                        }
                        if (g_match_msg_count == 3)
                        {
                            printf("The naming system does not appear to be consistent.\n No more match messages will be printed.\n");fflush(0);
                        }   
                        g_match_msg_count++;                     
                        break;
                    }
                }

            }
            if (result == 1) break;
        } 

    }

    free(query);  

    return result;
}


/* Scan the file, checking format and locating the clone, if any, which the 
        sequence matches in fpc, if it is a target sequence.
    Also build a tree of the lengths because blast doesn't provide them.  
    Also, compute the "sequence contig" number. This will be either

    A. The numeric order of the sequence in the file
    B. The number parsed from the .N or ContigN convention
    C. A mixture of both if convention B is not followed consistently.
*/
int collect_fasta_seq_data(char* path, FPC_TREE* ptree,struct bss_instance* pbsi, int target)
{
    FILE* f;
    int nseqs = 0;
    char line[1024];
    char query_id[1000];
    char query_orig[1000];
    SEQ_DATA* ptqd;
    int i;
    char* ptr, *ptr2;
    char seqname[100];
	char *tmpstr;
    int seqctg;
    int cidx = -1;
    gchar* m;

    int messresult;
    int munged;
    int gb_header;
    int order = 0;

    g_match_msg_count = 0;


    f = fopen(path,"r");
    if (f)
    {
        while (fgets(line,1024,f))
        {
            if (line[0] == '>')
            {
                order++;
                /* skip initial spaces and then parse to the first space */
                ptr = line;
                ptr++;
                while (isspace(*ptr)) ptr++;
                ptr2 = ptr;
                while (*ptr2 && !isspace(*ptr2) && *ptr2 != '\n')
                {
                    ptr2++;
                } 
                *ptr2 = 0;
                safe_strncpy(query_id,ptr,99);
                safe_strncpy(query_orig,ptr,199);

                if (strlen(query_id) >= BSS_FIELD_MAXSIZE)
                {
                    m = g_strdup_printf("Query name %s in file %s is too long (max length %d).", query_id,path,BSS_FIELD_MAXSIZE);
                    messout(m);
                    g_free(m);
                    return -1;
                }

                /* If it looks like a genbank header, tell the user we're going to pull
                    out the last field. We have to do this because the GB headers are too long. 
                    unless the program is megablast, in which case we have to match its 
                    idiotic parsing conventions */

                munged = 0;
                gb_header = 0;
                if (0 == strncmp(query_id,"gi|",3))
                {
                    gb_header = 1;
                    ptr = strrchr(query_id, '|');
                    if (pbsi->program != MEGABLAST)
                    {
                        m = g_strdup_printf("Found Genbank-style header %s. The last field (accession number) will be used as the query name.", query_id);
                        ptr++;
						tmpstr = strdup(ptr);
						strcpy(query_id,ptr);
						free(tmpstr);
						tmpstr = 0;

                    }
                    else if (target)
                    {
                        m = g_strdup_printf("Found Genbank-style header in target sequence %s. Megablast uses the third field (accn # with revision) as the sequence id.", query_id);
                        if ( (ptr = strrchr(query_id, '|')) )
                        {
                            *ptr = 0;
                        }
                        else
                        {
                            printf("Not a valid Genbank header:\n%s\nCannot continue.",line);
                            return -1;
                        }    
                        if ( (ptr = strrchr(query_id, '|')) )
                        {
                            ptr++;
                            memmove(query_id, ptr, strlen(ptr) + 1);
                        }
                        else
                        {
                            printf("Not a valid Genbank header:\n%s\nCannot continue.",line);
                            return -1;
                        }   
                    }
                    else
                    {
                        m = g_strdup_printf("Found Genbank-style header in query sequence %s. Megablast uses the first field (Genbank ID) as the sequence name.", query_id);
                        ptr = strchr(query_id, '|');
                        ptr++;
                        if ( (ptr = strchr(query_id, '|')) )
                        {
                            *ptr = 0;
                        }
                        else
                        {
                            printf("Not a valid Genbank header:\n%s\nCannot continue.",line);
                            return -1;
                        }
						tmpstr = strdup(ptr);
						strcpy(query_id,ptr);
						free(tmpstr);
						tmpstr = 0;
                    }
                    if (!g_gb_msg_shown && !Zbatch_flag)
                    {

                        messresult = messQuery(m);
                        g_free(m);
                        if (messresult == 0) 
                        {
                            return -1;
                        }
					}
                    g_gb_msg_shown = 1;

                }

                /* Set the seqctg number */
                if (!parse_seqctg(query_id,seqname,&seqctg))
                {
                    seqctg = order;
                    strcpy(seqname,query_id);
                }
                else
                {
                    if (seqctg < order)
                    {
                        seqctg = order;
                        if (0) //!g_seqctg_msg_shown && !Zbatch_flag)
                        {
                            m = g_strdup_printf("The sequence contig naming system .N or ContigN rmat >seqname.NN or >seqname.ContigN. Continue?");
                            messresult = messQuery(m);
                            g_free(m);
                            if (messresult == 0) 
                            {
                                return -1;
                            }
                        }
                        g_seqctg_msg_shown = 1;                          
                    }
                    else if (seqctg > order)
                    {
                        order = seqctg;
                    }
                } 
                /*see if the query id is too long for use as a marker */
                if (g_length_warn == 0 && strlen(query_id) > MARKER_SZ && !Zbatch_flag)
                {
                    m = g_strdup_printf("Query name %s is too long to be used as an FPC marker name (max %d). You will not be able to add its hits to FPC as markers. Continue?", query_id,MARKER_SZ);
                    messresult = messQuery(m);
                    g_free(m);
                    if (messresult == 0) 
                    {
                        return -1;
                    }  
                    g_length_warn = 1;                  
                }

                ptqd = fpc_tree_find_data(ptree,query_id);
                if (!ptqd)
                {
                    ptqd = (SEQ_DATA*)malloc(sizeof(SEQ_DATA));
                    memset(ptqd,0,sizeof(SEQ_DATA));
                    ptqd->cidx = -1;
                    strncpy(ptqd->name,query_id,99);
                    strncpy(ptqd->name_org,query_orig,199);
                    ptqd->munged = 1;
                    ptqd->seqctg = seqctg;
                    fpc_tree_insert_node(ptree, (void*)ptqd, query_id);
                }
                else
                {
                    if (!Zbatch_flag) messout("%s contains duplicated sequence name %s\n",path,query_id);
                    else printf("%s contains duplicated sequence name %s\n",path,query_id);fflush(0);
                    return -1;
                }
                /* if it's target side, try to match to a clone, unless we already matched the file */
                if (target)
                {
                    if (match_query_to_clone(seqname,&cidx,pbsi))
                    {
                        ptqd->cidx = cidx;
                        pbsi->seqs_matched++;
                    }
                }
                nseqs++;
            }
            else
            {
                assert(ptqd);
                for (i = 0; i < strlen(line); i++)
                {                   
                    if ( isalpha(line[i]))
                    {
                        ptqd->len++;
                    } 
                }  
            }
        }
        fclose(f);
    }
    return nseqs;
}

// the parameters are only used in batch mode
void do_search(parsetype type, int ctg, int from_end)
{
    struct bss_instance bsi;
    char temp[MAXPATHLEN];
    gchar* tempstr;
    FPC_TREE query_tree;
    FPC_TREE db_tree;
    int i;
    char* ptr;
    int n;
    struct stat s;
    gchar* m;
    DIR* dir;
    struct dirent* dent;
    int messresult;

    g_gb_msg_shown = 0;
    g_seqctg_msg_shown = 0;
    g_length_warn = 0;
    g_remove_bss_msg_shown = 0;

    memset(&bsi,0,sizeof(struct bss_instance));

    if (Zbatch_flag)
    {
        bsi.batch_search_type = type;
        bsi.batch_ctg = ctg;
        bsi.batch_fromend = from_end;
    }

    for (i = 0; i < BSS_NUMTYPES; i++)
    {
        bss_columns_checked[i] = 1;
    }

    bsi.filter_history = g_array_new(FALSE,FALSE,sizeof(BSS_FILTER));

    if (!dataloaded) 
    {
        fprintf(stderr,"No FPC file has been loaded.\n");
        return;
    }
    prepare_bss_results_dir();

    bsi.program = search_tool_type;
    bsi.split = split_bss_flag;

    ptr = dbDir;
    if (!strncmp(ptr,"./",2)) ptr += 2;
    strcpy(bsi.db_dir,ptr);

    ptr = qryDir;
    if (!strncmp(ptr,"./",2)) ptr += 2;
    strcpy(bsi.query_dir,ptr);

    ptr = qryFile;
    if (!strncmp(ptr,"./",2)) ptr += 2;
    strcpy(bsi.query_file,ptr);

    ptr = dbFile;
    if (!strncmp(ptr,"./",2)) ptr += 2;
    strcpy(bsi.db_file,ptr);

    safe_strncpy(bsi.query_suffix,qrySuffix,10);
    safe_strncpy(bsi.db_suffix,dbSuffix,10);
    if (!Zbatch_flag)
    {
        safe_strncpy(bsi.cutoff_str,(gchar*)gtk_entry_get_text(GTK_ENTRY(cutoff_entry)),10);
        safe_strncpy(bsi.xtraParm,(gchar*)gtk_entry_get_text(GTK_ENTRY(additional_entry)),99);
        strcpy(bsi.bss_subdir,new_bss_subdir());
    }
    else
    {
        strcpy(bsi.cutoff_str,cutoff_str);
        strcpy(bsi.xtraParm,xtraParm);
    }

    bsi.dbs = g_array_new(TRUE, FALSE, sizeof(gchar*));
    bsi.queries = g_array_new(TRUE, FALSE, sizeof(gchar*));

    bsi.old_filter_history = g_array_new(TRUE, FALSE, sizeof(char*));
    bsi.filter_history = g_array_new(TRUE, FALSE, sizeof(BSS_FILTER));

    if(!strcmp(bsi.query_file,"<all>"))
    {
        collect_filenames(bsi.query_dir,bsi.query_suffix,bsi.queries);
    }
    else
    {
        tempstr = g_strdup(bsi.query_file);
        g_array_append_val(bsi.queries,tempstr);
    }

    if(!strcmp(bsi.db_file,"<all>"))
    {
        collect_filenames(bsi.db_dir,bsi.db_suffix,bsi.dbs);
    }
    else
    {
        tempstr = g_strdup(bsi.db_file);
        g_array_append_val(bsi.dbs,tempstr);
    }

    if (bsi.dbs->len == 0)
    {
        fprintf(stderr,"No database files specified\n");fflush(0);    
        return;
    }
    if (bsi.queries->len == 0)
    {
        fprintf(stderr,"No query files specified\n");fflush(0);  
        return;  
    }    

    /* scan all the target files, check for unique names, and build a table of seq lengths */

    init_fpc_tree(&db_tree);

    bsi.num_initial_targets = 0;
    for (i = 0; i < bsi.dbs->len; i++)
    {
        printf("\n");fflush(0);
        sprintf(temp,"%s/%s",bsi.db_dir,g_array_index(bsi.dbs,char*,i));
        printf("Computing prefix of clones in: %s\n",temp);fflush(0);
        bsi.seqs_matched = 0;
        n = collect_fasta_seq_data(temp,&db_tree,&bsi,1);
        if (-1 == n)
        {
            clear_fpc_tree(&db_tree,delete_seq_data);
            return;
        }
        printf("%d sequences (out of %d) could be matched to a clone\n",bsi.seqs_matched,n);fflush(0);

        if (!Zbatch_flag && ((float)bsi.seqs_matched)/((float)n) < .5)
        {
            m = g_strdup_printf("The names of sequences in %s do not match well with clone names in FPC. Continue?", temp);
            messresult = messQuery(m);
            g_free(m);
            if (messresult == 0) 
            {
                clear_fpc_tree(&db_tree,delete_seq_data);
                return;
            }   
        }

       
        bsi.num_initial_targets += n;
    }

    for (i = 0; i < bsi.queries->len; i++)
    {

        sprintf(bsi.query_file,"%s",g_array_index(bsi.queries,char*,i));

        sprintf(temp,"%s/%s",GETCURDIR,"BSS_results");
        if (strlen(bsi.bss_subdir) > 0)
        {
            sprintf(temp,"%s/BSS_results/%s",GETCURDIR,bsi.bss_subdir);   
        }
        if (!Zbatch_flag)
        {
            sprintf(bsi.bssfile,"%s/%s",temp,bsi.query_file);
            if ( strrchr(bsi.query_file,'.') && (ptr = strrchr(bsi.bssfile,'.'))  )
            {
                *ptr = 0;
            }
            if (bsi.dbs->len == 1)
            {
                strcat(bsi.bssfile,".");
                strcat(bsi.bssfile,g_array_index(bsi.dbs,char*,0));
                if ( strrchr(g_array_index(bsi.dbs,char*,0),'.') && (ptr = strrchr(bsi.bssfile,'.')) )
                {
                    *ptr = 0;
                }
            }
            else
            {
                strcat(bsi.bssfile,".");
		        ptr = bsi.db_dir;
		        if (strrchr(bsi.db_dir,'/'))
		        {
		           ptr = strrchr(bsi.db_dir,'/');
                   ptr++;
		        }
                strcat(bsi.bssfile,ptr);
                if (strlen(bsi.db_suffix) > 0)
                {
                    strcat(bsi.bssfile,".");
                    strcat(bsi.bssfile,bsi.db_suffix);
                }
            }
            strcat(bsi.bssfile,".bss");
        }
        else
        {
            strcpy(bsi.bssfile,batch_bss_file);    
        }

        if (stat(bsi.bssfile, &s) == 0 && !Zbatch_flag && !g_remove_bss_msg_shown)
        {
            m = g_strdup_printf("%s already exists. Previously-existing BSS results files will be removed. Continue?", bsi.bssfile);
            messresult = messQuery(m);
            g_free(m);
            g_remove_bss_msg_shown = 1;
            if (messresult == 0) continue;
            if (S_ISDIR(s.st_mode))
            {
               dir = opendir(bsi.bssfile);
               if (dir != NULL) 
               {
                  while((dent = readdir(dir)) != NULL)
                  {                     
                      if (strcmp(".",dent->d_name) && strcmp("..",dent->d_name))
                      {
                         sprintf(temp,"%s/%s",bsi.bssfile,dent->d_name); 
                         remove(temp);
                      }
                   }
                }
                closedir(dir);
                rmdir(bsi.bssfile);               
            }
            else
            {
                remove(bsi.bssfile);
            }
            if (stat(bsi.bssfile, &s) == 0)
            {
                printf("Unable to remove %s; remove it or choose a different name to continue.\n",bsi.bssfile);fflush(0);
                return;
            }
        }
    

        sprintf(temp,"%s/%s",bsi.query_dir,bsi.query_file);
        printf("Collect sequence lengths from file:%s\n",temp);fflush(0);
        init_fpc_tree(&query_tree);
        bsi.num_initial_queries = collect_fasta_seq_data(temp,&query_tree,&bsi,0);
        if (-1 == bsi.num_initial_queries)
        {
            clear_fpc_tree(&query_tree,delete_seq_data);
            return;
        }
        switch  (bsi.program)
        {
            case BLAST:
            case MEGABLAST:
                blast_search(&bsi,&query_tree,&db_tree);
                break;
            case BLAT:
                blat_search(&bsi,&query_tree,&db_tree);
                break;
            default:
                assert(0);
        }
        clear_fpc_tree(&query_tree,delete_seq_data);
    }

    clear_fpc_tree(&db_tree,delete_seq_data);

    printf("Number of files processed: %d; Total retained hits:%d\n",bsi.queries->len, bsi.retained_hits);
    printf("On average, each file had %d hits.\n", bsi.total_hits/bsi.queries->len);
}


void blat_search(struct bss_instance* pbsi, FPC_TREE* pqtree, FPC_TREE* pdbtree)
{
    char listfile[MAXPATHLEN];
    char outFile[MAXPATHLEN];
    char cmd_str[1024];
    int i;

    FILE* lfp;

    /* make the database list file  */

    sprintf(listfile,"%s/blat.list",pbsi->db_dir);
    lfp = fopen(listfile,"w+");
    if (!lfp)
    {
        fprintf(stderr,"Unable to create %s\n",listfile);fflush(0);
        return;
    }
    for (i = 0; i < pbsi->dbs->len; i++)
    {
        fprintf(lfp,"%s/%s\n",pbsi->db_dir,g_array_index(pbsi->dbs,char*,i));
    }
    fclose(lfp);

    sprintf(outFile,"%s.blat",pbsi->bssfile);
    pbsi->blat_score = strtod(pbsi->cutoff_str,0);
    sprintf(cmd_str,"blat %s %s/%s %s -minScore=%d %s", listfile,pbsi->query_dir,pbsi->query_file,pbsi->xtraParm,pbsi->blat_score,outFile);

    fpc_exec(cmd_str);
    sprintf(cmd_str,"rm %s",listfile);
    system(cmd_str);
    parse_blat_output(pbsi,outFile,pqtree,pdbtree);
    reset_bss_results();
}
void blast_search(struct bss_instance* pbsi, FPC_TREE* pqtree, FPC_TREE* pdbtree)
{
    char outFile[MAXPATHLEN];
    char nalFile[MAXPATHLEN];
    FILE* nalfp;
    int i;
    char temp[MAXPATHLEN];
    char seq_path[MAXPATHLEN];
    char cmd_str[1024];
    struct stat s;
    char formatdb_dir[MAXPATHLEN];
    char *seqfile;
    time_t seq_mtime;

    /* Make the formatdb directory if it doesn't exist */

    sprintf(formatdb_dir,"%s/formatdb",pbsi->db_dir);
    if (0 != stat(formatdb_dir,&s))
    {
        sprintf(cmd_str,"mkdir %s",formatdb_dir);
        fpc_exec(cmd_str); 
    }
    /* Format the db files */

    printf("Formatting BLAST databases\n"); fflush(0);
    for (i = 0; i < pbsi->dbs->len; i++)
    {
        /* get the mod time of the sequence file */
        seqfile = g_array_index(pbsi->dbs,char*,i);
        sprintf(seq_path,"%s/%s",pbsi->db_dir,seqfile);
        assert(0 == stat(seq_path,&s));
        seq_mtime = s.st_mtim.tv_sec;

        sprintf(temp,"%s/%s.nsq",formatdb_dir,seqfile);

        if (0 != stat(temp,&s) || s.st_mtim.tv_sec < seq_mtime)
        {
            printf("%s needs to be formatted\n",seqfile);fflush(0);
            sprintf(cmd_str,"formatdb -i %s -p F  -n %s/%s",seq_path,formatdb_dir,seqfile);
            printf("%s\n",cmd_str);fflush(0);
            fpc_exec(cmd_str);
            strcat(temp,".nsq");        
            if (0 == stat(temp,&s))
            {
                printf("Failed to create %s; aborting.\n",temp);fflush(0);
                return;
            }    
        }    
    }

    /* make the .nal file for the group */

    sprintf(nalFile,"%s/bss.nal",formatdb_dir);
    nalfp = fopen(nalFile,"w+");
    if (!nalfp)
    {
        fprintf(stderr,"Unable to create %s\n",nalFile);fflush(0);
        return;
    }
    fprintf(nalfp,"TITLE: BSS Search\n");
    fprintf(nalfp,"DBLIST ");
    for (i = 0; i < pbsi->dbs->len; i++)
    {
        fprintf(nalfp," %s",g_array_index(pbsi->dbs,char*,i));
    }
    fprintf(nalfp,"\n");
    fclose(nalfp);



    if (pbsi->program == MEGABLAST)
    {
        sprintf(outFile,"%s.megablast",pbsi->bssfile);
        sprintf(cmd_str,"megablast -V F -d %s/bss -i %s/%s -o %s -e %s %s -D 3", 
            formatdb_dir, pbsi->query_dir,pbsi->query_file, outFile, pbsi->cutoff_str,  pbsi->xtraParm);
    }
    else
    {
        sprintf(outFile,"%s.blast",pbsi->bssfile);
        sprintf(cmd_str, "blastall -p blastn -d %s/bss -i %s/%s -o %s -e %s %s -m 8",
                     formatdb_dir, pbsi->query_dir,pbsi->query_file, outFile, pbsi->cutoff_str, pbsi->xtraParm);
    }
    fpc_exec(cmd_str);
    parse_blast_output(pbsi,outFile,pqtree,pdbtree);
    reset_bss_results();
}

// prepare_bss_results_dir
//
// Test for BSS_results directory, and create it if necessary.
//
static void
prepare_bss_results_dir(void)
{
   int here;
   struct stat stat_buf;
   const gchar* subdir;

   here = open(".", O_RDONLY);
   g_assert(here >= 0);
   chdir(GETCURDIR);
   if (stat("BSS_results", &stat_buf) != 0)
       mkdir("BSS_results", 0777);
   else if (!S_ISDIR(stat_buf.st_mode)) {
       unlink("BSS_results");
       mkdir("BSS_results", 0777);
   }
   subdir = new_bss_subdir();
   if (subdir[0] != '\0') {
       chdir("BSS_results");
       if (stat(subdir, &stat_buf) != 0)
           mkdir(subdir, 0777);
       else if (!S_ISDIR(stat_buf.st_mode)) {
           unlink(subdir);
           mkdir(subdir, 0777);
       }
   }
   fchdir(here);
}

